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ABSTRACT 

The shallow water equations are a set of coupled hyperbolic partial differential equa- 
tions. As analytical solutions for these equations are not possible for practical cases 
of open channel flow, they are solved using numerical techniques. Finite difference 
and finite element methods have been applied in the past for solving most of the 
practical cases. However, Finite Volume Methods are becoming popular these days 
because they can be applied to any arbitrary physical domain (like finite element 
methods) and are easy to implement (unlike finite element methods). In this study, 
an overlapped control volume (OCV) method is developed for solving the two-, 
dimensional shallow water flow equations. The OCV method is verified by solving 
steady two-dimensional supercritical flow problems. The test cases solved are chan- 
nel transitions, the results of which are compared with the solutions obtained by 


employing a finite difference method. 
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NOTATION 


CN Courant Number 

g acceleration due to gravity, m/s 2 

F Froude number 

h flow depth , m 

i subscript for space node in X direction 

j subscript for space node in Y direction 

n subscript for time level 

S ox channel bottom slope in x-direction 

S oy channel bottom slope in y-direction 

Sf x friction slope in x-direction 

Sf y friction slope in y-direction 

n Manning’s roughness coefficient 

t time 

u flow velocity in x-direction 

v flow velocity in y-direction 

V magnitude of resultant velocity 

x . space coordinate 
y space coordinate 



Greek letters 


a x angle between channel and the ^-direction 

ot y angle between channel and the y-direction 

k dissipation constant used in the numerical scheme 

At time interval 

Ax distance increment in x— direction 

A yi distance increment in y— direction at x — i 

0 angle between the wall and the x- axis 

S inclination of velocity vector at interior node adjacent to boundary 

1 / scaling parameter used in artificial viscosity procedure 

u> weighting coefficient 
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Chapter 1 


Introduction 


Shallow water flow equations are generally used for simulating most of the un- 
steady open channel flows including flood waves, dambreak flows, flows in estuaries, 
etc. They are also used often to represent two-dimensional supercritical flows with 
oblique waves when the shallowness parameter is very small. These equations are 
a set of coupled hyperbolic partial differential equations . Analytical solutions of 
these equations are available only for idealized cases; therefore, they are numerically 
solved for most of the practical applications. 

Many numerical methods are available for the solution of the shallow water flow 
equations. These numerical techniques can be broadly classified into the following 
types: 

1. Finite-difference Methods, 

2. Finite-element Methods and 

3. Finite-volume Methods. 

Historically, Finite-difference methods are the oldest and they are popular even 
today. Finite-difference methods such as the Priessman scheme, Lax-Wendroff 
scheme etc. have been traditionally used for solving the one-dimensional shal- 
low water flow equations (Chaudhry, 1993). Fennema and Chaudhry (1990) and 
Garcia and Kawahita (1986) introduced the MacGormack scheme for solving the 
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two-dimensional shallow water flow equations. Since then, this method has been 
extensively used for solving various two-dimensional open-channel flow problems 
(Jimenez and Chaudhry 1988 , Dammuller et al. 1989 , Bhallamudi and Chaudhry 
1992). Fennema and Chaudhry (1990) have also introduced the implicit Beam and 
Warming method for solving two-dimensional dambreak flow problems. Casulli 
(1990) developed ADI, semi-implicit and fully implicit splitting Finite-difference 
methods for solving the two-dimensional shallow water flow equations. Classical 
second-order accurate finite-difference methods result in dispersive errors when ap- 
plied to flows with shocks. These dispersive errors manifest in terms of higher-order 
oscillations near the shock front. Artificial viscosity methods are usually applied to 
smoothen these higher-order oscillations (Chaudhry, 1993). Recently, Total Varia- 
tion Diminishing (TVD) schemes (Garcia-Navvarro et al. 1992, Alcrudo et al. 1992) 
and Essentially Non-Oscillatory (ENO) schemes (Yang et. al. 1993 , Nujic 1995) 
were introduced for solving the shallow water flow equations. These techniques have 
a sounder theoretical basis than the artificial viscosity procedures for smoothing the 
higher-order oscillations. Savic and Holly (1993) developed the modified Godunov 
method for solving a one -dimensional dambreak Hood wave. This method is very 
useful for solving mixed flow regimes and for flows with strong shocks, without the 
introduction of artificial viscosity; however, it has not been used for two-dimensional 
flows. 

The main advantage of any finite-difference method is the ease of implementa- 
tion. They are best suited for one-dimensional applications. However, their applica- 
tion to two-dimensional flows requires approximation at the boundaries if the flow 
domain is not rectangular. Alternatively, Coordinate transformation techniques may 
be used for converting the non-rectangular physical domain into a rectangular com- 
putational domain (Dammuller et al. 1989, Bhallamudi and Chaudhry 1992, Bellos 
et al. 1991). Bhallamudi and Chaudhry (1992) used a global transformation of the 
governing equations while Bellos et al. (1991) used the transformation of equations 
locally for the non-orthogonal finite-difference grid. This requirement of coordinate 
transformation makes the application of finite-difference methods difficult in case 
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of two-dimensional open channel flows in general non-rectangular domains. 

Finite-element methods for shallow water flow equations have been developed es- 
sentially to overcome the limitations of finite-difference methods for two-dimensional 
flow in arbitrary domains. Katapodes (1984) developed a dissipative Galerkin 
scheme for simulating two-dimensional surges and shocks in open channels. The 
parasitic waves in the vicinity of the discontinuity, commonly present in the clas- 
sical Galerkin solution, are selectively dissipated in this technique. The method is 
based on discontinuous weighting functions which introduce upwind effects in the 
solution. Akanbi and Katapodes (1988) modified the above technique for appli- 
cation to a deforming coordinate system to simulate flood waves propagating on 
dry bed. The U.S.Army Corps of Engineers (1985) developed a commercial finite- 
element package for the simulation of unsteady flow in estuaries and bays. This 
model incorporates a constant eddy viscosity model for turbulence, which helps in 
dissipating the higher-order oscillations (Chaudhry et al. 1988). Lynch and Gray 
(1978) and Kawahara and Umetsu (1986) developed finite-element methods for mov- 
ing boundary problems in river flows. Herrling (1982) developed a finite-element i 
model based on coupling of one- and two -dimensional finite elements for simulat- 
ing moving boundary problems in estuaries. Leclere et al. (1990) also developed 
a finite-element model for shallow water flow in estuaries with moving boundaries. 

This model is based on a fixed spatial mesh with elements becoming dry or wet j 

in a continuous manner. Recently, Galland et al. (1991) developed a commercial j 

software package called TELEMAC for simulation of shallow water flows using a j 
finite-element method. Although finite-element methods can be applied to any ir- j 
regular flow domains, algorithms based on these techniques are quite complex and | 

require much more computational power than the finite-difference methods for their j 

application to even simple problems. j 

Finite-volume techniques offer a viable alternative to the finite-element meth- j 

ods for solving the fluid flow problems (Peyret and Taylor 1983). They com- j 

bine the flexibility of handling the complex geometries, intrinsic to FEM , with j 

the simplicity of FDM. Putti et al. (1990) used a triangular finite-volume tech- j 
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nique for solving the ground water solute transport equation. Verma and Eswaran 
(1996) developed an Overlapping Control Volume (OCV) method for the solution 
of two-dimensional Convective-Diffusion equation. Later, Verma et al. (1996) ex- 
tended this method for simulating two-dimensional solute transport in ground water. 
Alcrudo and Garcia-Navvarro (1993) introduced a high resolution Godunov-type 
scheme in finite-volumes for the two-dimensional shallow water flow equations. The 
method proposed by Alcrudo and Garcia-Navvarro (1993) can be applied to any ar- 
bitrary flow domain without the necessity of coordinate transformation. 

In this study, an attempt is made to extend the OCV technique developed by 
Verma and Eswaran (1996) for application to the shallow water flow equations. 
Unlike the Convective-Diffusive equation solved by Verma and Eswaran (1996) , the 
shallow water flow equations do not contain the diffusive terms. On the other hand, 
the complications arise due to the coupled nature of the Continuity and Momentum 
equations. The proposed method is verified by solving two-dimensional supercritical 
flow problems in channel transitions. The numerical results obtained obtained using 
the OCV technique compare satisfactorily with earlier results obtained using the 
Finite-difference method. 

The governing equations for the shallow water flow are presented in Chapter 2 
followed by a description of the OCV technique in Chapter 3. The OCV technique 
is validated using earlier numerical results obtained by a finite-difference method 
for channel transitions. Concluding remarks are presented in chapter 5. 



Chapter 2 


Governing Equations 

The two-dimensional shallow water flow equations represent the laws of conservation 
of mass and momentum for flow in open-channels. These equations are derived from 
the three-dimensional equations of motion by integrating them over the flow depth. 
The shallow water equations in cartesian system are written as, 

Ut + F s + G y + S= 0 ( 2 . 1 ) 

where 

h 

hu 
hv 

uh 

F = u 2 h + \<jh l 

i uvh 
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G = 


( 


vh 

uvh 

v 2 h + \gh '? 


( 


S = 


0 

-gh(S ox - Sf X ) 

Soy Sfy) 




subscripts t , x and y in Eqn 2.1 represent the partial derivatives with respect to 
time, longitudinal direction and lateral direction, respectively. 

In the above equation, 
t = time 

u = velocity along the longitudinal direction 

v = velocity along the lateral direction 

h = water depth 

g = acceleration due to gravity 

S ox = channel slope along the longitudinal direction 

S oy = channel slope along the lateral direction 

S fx = friction slope along the longitudinal direction 

Sfy = friction slope along the lateral direction 


The bottom slopes are estimated as: 


S ox = sin a x 

(2.2) 

S oy = sin a y 

(2.3) 


a x ,a y are the angles between the bottom of the channel and the longitudinal and 
lateral direction, respectively. 


j 


1 
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Friction slopes may be estimated using the Manning’s or Chezy’s equation. 


5 fx 


Sfy = 


+ v 2 

w* 

n 2 Vy/u 2 + V 2 


h 4 / 3 


(2.4) 

(2.5) 


where n = Manning’s roughness coefficient. 


2.1 Assumptions 

Shallow water flow equations are based on the following assumptions : 

1. The pressure distribution in the vertical direction is assumed to be hydrostatic 
i.e. the acceleration in the vertical direction is assumed to be negligible. 

2. Depth averaged values are used in two-dimensional modelling. 

3. Fluid is assumed to be incompressible and the density is same throughout. 

4. Velocity distribution is uniform along the depth. 

5. The channel bottom is rigid and the bottom slope is small. 

6. Only shear stresses due to horizontal velocity components are significant, other 
turbulence effects are neglected. 

7. Frictional resistance is approximated by the Chezy s equation or Manning s 
equation. 

Although shallow water equations are applicable for most of the open-channel flows 
under subcritical conditions, their application to supercritical flows is limited be- 
cause of the assumptions 1 and 2. Therefore, these assumptions are discussed briefly 

in the following paragraphs. 



8 


2.1.1 Hydrostatic Pressure distribution: 

The assumption of hydrostatic pressure distribution is correct if the streamlines are 
parallel and straight. As long as high curvatures are not present, the pressure distri- 
bution is almost hydrostatic. Ligget (1975) had shown that the assumption is valid 
as long as a shallowness parameter Hq/Iq (Jiq is water depth and Iq is characteristic 
length) is small. Later, Jimenez (1987) concluded that the shallow water theory 
reasonably represents the supercritical flow if depth to width ratio is of the order of 
0.1 and Fronde number is not close to unity. The error is of the order of 20% if the 
depth to width ratio is increased to 0.2. This error is manifested in the wavelength 
of the resulting wave pattern. The hydrostatic assumption is not valid in the vicin- 
ity of the shock and some flow details are lost at the shock location. However the 
overall results are adequate for engineering purposes (Cunge 1975). 

2.1.2 Turbulent effects and Depth averaging: 

The depth integration of the three-dimensional equations of motion produce higher 
order terms known as effective stresses which are not included in the Eq(2.1). Effec- 
tive stresses constitute laminar viscous stresses, turbulent stresses and stresses due 
to depth averaging of advective terms. Laminar viscous stresses may be neglected 
but the other two are important. Consideration of turbulent stresses requires the 
use of a turbulence closure model which expresses stresses as a function of the main 
flow variables. Rastogi and Rodi (1978), Puri and Kao (1984) and Leschziner and 
Rodi (1979) used K - e models for describing turbulence. Vreugdenhil and Wi- 
jbenga (1982) used the constant eddy viscosity concept. However, the above closure 
models have been tested only for subcritical flows and no information is available for 
supercritical flows. Inclusion of turbulence effects is beyond the scope of the present 
study. 

Ippen and Harleman (1950) studied the effect of non-uniform velocity distribu- 
tion and concluded that the effect is negligible.The measurement of velocity across 
the depth on the front and back side of the jump indicated that the momentum cor- 
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rection factor ranged from 1.007 at F r = 6.3 to 1.015 at F r = 2.0. It should be noted 
that though the effect of magnitude variation of velocity may not be significant, yet 
the variation of velocity across the depth may be significant (Kumar, 1996). 

In the present study, it is assumed that the shallow water equations are appli- 
cable for supercritical flows also. Their numerical solution is presented in the next 
chapter. 



Chapter 3 


Overlapping Control Volume 
(OCV) Method 


The shallow water flow equations are nonlinear hyperbolic partial differential equa- 
tions and their closed form solutions are available only for very simplified cases. 
Therefore, they are numerically integrated using a control volume approach. An 
Overlapping Control Volume (OCV) method adapted from Verma and Eswaran 
(1996) and Verma et al. (1996) is used for this purpose. However, the flux cal- 
culation at the control surface of the control volume is different in the present ap- 
plication. Also, the partial discretization for time is done using a MacCormack 
like scheme (MacCormack, 1969). The proposed scheme is an explicit, two-step 
predictor-corrector scheme which is second order accurate, in time and space, and 
is capable of capturing shocks without isolating them. The method is applicable on 
non-orthogonal grids. 


3.1 Finite- Volume Formulation 

The physical domain over which the solution is sought is discretized into a structured 
non— orthogonal grid as shown in Fig 3.1a . A Control Volume as shown in Fig 3.1b is 
considered. The Control Volume is labelled by the index of the control node i.e., as 
(i j). This type of Control Volume uses the grid point coordinates to form the Control 
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Volumes. Therefore, unlike in the usual Control Volume approach, it does not 
involve the determination for any of the intermediate points. However, the adjacent 
Control Volumes will overlap each other: hence the name “Overlapping Control 
Volume technique. The governing equation for the control volume is obtained by 
integrating the Eq.(2.1) over the control volume and applying the Gauss divergence 
theorem. 


J c J^r dA + j> c {Fn x + Gn y )dl = J J SdA (3.1) 

where, dl is an elemental length on the boundary (Control Surface) of the Control 
Volume, n x and n y are the direction cosines of the local outward unit vector on the 
boundary in x and y directions, respectively, and dA is the elemental area of the 
Control Volume. 


3.2 Solution of Governing Equation 

Equation (3.1) is solved for any control volume (i,j) using an explicit predictor- 
corrector approach. Discretization of Eq (3.1) for an explicit solution for time level 
t + At is explained in the following sections. 

3.2.1 Predictor Step 

Equation (3.1) is partially discretized in the predictor step as, 

77 * _ TJ n 

= S?jA s - (Flux) n (3.2) 

where 

(Flux) n = I {Fn x + Gn y )dl (3.3) 

As is the area of the Control Volume and A t is the computational time step. 
Subscript (i,j) indicates the Control Volume and the superscripts * and n in Eqs. 
(3.2) and (3.3) denote the evaluation of the terms at the predictor level and the time 
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level t , respectively. The area of the control volume (.A s ) is calculated using the. 
formula 

A* ~ 0 - 5 [(*W+i ~ ~ Vi+ 1 , 1 ) ~ ~ x i+ ij)(yij+i ~ Vi+i,j) 

+ {xi,j-i ~ x i-\j){yi,j+i - yi-ij) - (x iij+l - - Vi-hj)] 

where etc., are the coordinates of the neighbouring points shown 

in Fig. 3.1b. The contour integration for the term Flux in Eq.(3.2) is counter- 
clockwise. Discretization of this term is explained in detail below. 


3.2.2 Continuity Equation: 

The U variable in Eq.(3.2) is h for the continuity equation. For determining the 
predicted flow depth, h*, Flux is discretized as 


Flux = J2 (u (fc) Ay (fc) -u (fc) Ax w ) 

k=l 

= (3.4) 

k= 1 

where the superscript ( k ) refers to the edges of the control volume(see Fig 3.1b ). 
For the edge ( k ), k= 1,2 and 3, the approximation used is 


u (fc) = 0.5(u* + u k+ i) (3.5) 

v {k) = 0.5(u fc + u fc+1 ) . (3.6) 

A = (Vk+i ~ Vk) ( 3 - 7 ) 

Ax (fc) = (xfe+i - Xk) ( 3 - 8 ) 


where the subscript k refers to the grid point number as shown in Fig 3.1 b. For 
k= 4, u k+h x k+ i etc. are replaced by u u x i respectively, in the above equations. 

To incorporate upwinding, h « in Eq.(3.4) is approximated at the mid-point 
of control-surface k by interpolation within the upwind control volume adjacent 
to the surface. For example, referring to Fig. 3.1 b , if the flow is entering the 
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control volume (i,j) across face 1 ,i.e., J (1) is negative, then h W is approximated by- 
interpolation within the control volume for node — — The values at the grid 

points constituting the control volume (z - 1, j - 1) are used for determining h 
On the other hand, if is positive, the values at the grid points of the control 
volume (z, j) are used to obtain h^. The interpolation scheme to obtain the 
value at the face k is explained below. 

Consider face 3 of the control volume (z, j) shown in Fig 3.2. h ^ at the mid- 
point of face 3 is obtained using the grid-point values of h for the control volume 
(i,j) if is positive. For this purpose, Verma and Eswaran (1996) and Verma 
et al. (1996) suggested an interpolation scheme based on finite-element type shape 
functions. However, it was found during the code development that such an inter- 
polation scheme did not give good results in case of shallow water equations. This 
was particularly so when a MacCormack type discretization was used for the time 
stepping. Therefore the following procedure was adopted in the present study. 



the diagonal joining the grid points 4 and 5 (see Fig 3.2), then /z (3) is obtained as 

= 0.5 (/i/i + /is) (3-9) 
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If the velocity vector at face 3 cuts the diagonal joining the grid points 3 and 5, 
then is obtained as 

h (3) = 0.5 (/i 3 + h$) (3.10) 

Similar procedure is adopted while determining h ^ and Substitution 

of Eqs (3.4)-(3.10) in Eq.(3.2) results in an explicit equation for the solution of 
predicted flow depth at ( i,j ) i.e., h*p since all other terms are evaluated at time 
level ”t” at which the solution is known. 


3.2.3 Momentum Equation 

The U variable in Eq. 3.2 is ( uh ) for solving the momentum equation in x-direction. 
For determining the predicted value of uh* , the Flux term is discretized as 

Flux = 

fc=l k = 1 1 

= + (3-ii) 

fc=l k = 1 1 

in which, the definition of superscripts is same as before. The approximation 
for u( k \ y( k \ Ax( k \ A .yW and h ^ are also same as described for the continuity 
equation. The term h W is obtained as 

/i' (fc) = 0.5 (h k + h k+l ) (3-12) 

for k= 1,2 and 3. is obtained by replacing Ajb-t-i by h x in Eq.(3.12). The 
upwinding and the interpolation scheme to obtain ( uh at the mid— point of the 
face k are similar to the procedures described for obtaining in the continuity 
equation. Substitution of Eqs. (3.5)-(3.12) in Eq.(3.2) results in an explicit equation 
for the solution of {uh)*j since all other terms are evaluated at the time level t. The 
predicted value of x-velocity at (z, j ) i.e., u*j can then be detennined as, 
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<i = 


\ uh )ii 


/>*. 

n hj 


(3.13) 


The procedure for determining v*j using the momentum equation in y-direction 
is similar to the procedure described for determining u*j. The Flux term in this 
case is discretized as 


Flux = AxW] 

k=l k = 1 2 

= (3.14) 

k = 1 *=1 2 

3.2.4 Corrector step 

The predicted values of depth and velocities i.e. /i*, tt* and v* are used in a corrector 
step to obtain the corrected values of depth and velocities. Superscript ** in the 
following equations represents the values of h, u and v at the end of the corrector 
step. Eq. (3.1) is partially discretized in the corrector step as 

[/*.* - c/» 

— M A S = S*jA s - (Flux)* (3.15) 

The discretization of the term Flux* in Eq.(3.15) is similar to the discretization 
of the term Flu x n in Eq.(3.2) except that downwinding is incorporated instead of 
the upwinding and predicted values are used instead of the values at time level 
t. This means that h^ k \ (uh)^ and ( vh )M in Eqs. (3.4), (3.11) and (3.14) are 
approximated at the mid-point of the control-surface k by interpolation within 
the downwind control volume adjacent to the surface. For example, referring to 
Fig. 3.1b, if the flow is entering the control volume (i,j) across face 1 , i.e. is 
negative, then h (1) , (uh) (1) and (vh) (1) are approximated by interpolation within the 
control volume for node (i,j)- The predicted values at the grid points constituting 
the control volume (i,j) are used for determining h (1) , (uh) (1) and (u/i) (1) . Rest 
of the procedure to determine h**, u** and v** is same as that described for the 
predictor step in the previous section. 
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3.2.5 Solution at the time level t + At 

The values of flow depth and velocity at the time level t + At are obtained by taking 
an average of the predicted and the corrected values of U. 


= 0-5 (U!j + u;j) (3.16) 

3.3 Initial and Boundary Conditions 

To start the unsteady state computations, the values of three primary variables h , u 
and v at time = 0 are specified at all the grid points. In the present application of 
the proposed numerical scheme, unsteady flow equations are used in a false transient 
approach to determine the steady state solutions. In the false transient approach, 
uniform velocity and uniform flow depth equal to the upstream boundary values 
are specified at all the nodes as the initial condition. It should be noted that any 
arbitrary values may be assumed as the initial condition since only the final steady 
state solution is of interest. 

The determination of u, v and h during the unsteady computations requires the 
application of the boundary conditions at the boundaries of the flow domain. Hyper- 
bolic equations are sensitive to the numerical treatment of the boundary conditions 
because errors introduced at the boundaries propagate and reflect throughout the 
domain and in many cases the instability may result (Anderson et. al. 1984). Three 
type of boundaries are encountered in the present study: 

1. The open (flow) boundaries. 

2. The symmetry boundary. 


3. Solid side wall boundary. 
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3.3.1 Flow Boundaries 

Inflow and outflow boundaries are open boundaries where flow can enter or leave 
the computational domain. The specification of boundary conditions depends on 
whether the flow is subcritical or supercritical. ( Stoker, Verbroom et al.). In 
the present study, the proposed numerical scheme is applied for solving the two- 
dimensional supercritical flow. Therefore, three boundary conditions are to be spec- 
ified at the inflow boundary and none at the outflow boundary. In the present study 

u, v and h are specified at the inflow boundary of the channel. These quantities 
at the downstream nodes are obtained by extrapolation from the already computed 
values at the interior nodes. Since the flow is supercritical, the extrapolation affects 
only a few points at the downstream end and the errors will not propagate upstream. 

3.3.2 symmetry boundary 

Reflection procedure is implemented at a symmetry boundary (Roache 1972). • All 
non-conservative flow variables other than the normal velocity ( u and h) are specified 
as even functions with respect to the symmetry line while the normal velocity is 
specified as an odd function so that the average normal velocity at the symmetry 
boundary is zero. The reflection technique is exact for a straight symmetry line. 

3.3.3 Solid Side wall boundary 

Since all the stresses other than the bottom stresses are neglected, the slip condition 
is the proper boundary condition for the solid boundaries. The basic requirement 
at a solid wall is that no flow should take place normal to the boundary, which is 
expressed as 

tan0 = - (3.17) 

u 

where 6 is the single between the wall and the x—SLxis. Therefore , the resultant 
velocity at a solid wall is tangent to it. Abbott (1971), Roache (1972) and Anderson et 

al. (1984) discussed several wall boundary techniques. Jimenez(1987) tested several 
techniques for the case of steady supercritical flow. The reflection procedure is used 
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in the present study. It should be noted that this procedure for the solid side walls is 
only approximate. Implementation of the reflection procedure depends on the sign 
and the magnitude of the values of 9 and S (S = inclination of the velocity vector 
at the interior node adjacent to the boundary). 

One of the possible scenarios for the case of channel expansion is explained here 
for illustration. Similar results can be obtained for other cases. Referring to Fig. 
(3.3) , the flow depth and the magnitude of resultant velocity at the imaginary 
reflection point are same as that of their values at the corresponding interior grid 
point. The direction of the flow velocity is determined such that the normal velocity 
at the wall is zero (Bhallamudi and Chaudhry,1992). 

hfic = hint (3.18) 

Vfic = Vint (3.19) 

where subeript fic refers to the values at the ” fictitious” reflection point and 
subscript int refers to the values at the interior point adjacent to the wall. If 9 is 
the angle between the wall and the z-axis and S is the angle between the resultant 
velocity at the interior point and the x- axis ,then, the velocity components uj ic 
and v^c at the reflection point are given by , 


Vint \] 'U’int ^ int 

(3.20) 

Ufic = Vint COS (20 - 6) 

(3.21) 

Vfic = Vint sin(20 - 5) 

(3.22) 


3.4 Stability Condition 

The proposed numerical scheme uses an explicit time stepping procedure. Therefore, 
Courant-Friedrichs-Lewy(CFL) condition has to be satisfied for numerical stability. 
Shallow water equations are nonlinear and hence this condition, though not strictly 
applicable, is used as a guideline. The condition for two dimensional flows is given 

by (Roache 1972), 



19 


A t < 


CN 






+ 


(3.23) 


A3/ 2 


where V is the velocity at the grid point and CN is the courant number. CN should 
be less than one for obtaining numerically stable results. 


3.5 Artificial Viscosity 

The proposed scheme is second order accurate in both space and time and will result 
in dispersive errors. These dispersive errors cause high-frequency oscillations near 
steep gradients. To smoothen these oscillations, a procedure developed by Jameson 
et al. (1981) is employed. This procedure smoothens regions of large gradients while 
leaving smooth areas relatively undisturbed. The values of the variables at the new 
time computed by the proposed method are modified using the following algorithm. 
A scaling parameter u is first determined based on the water level variation. 


_ hi- ij)| 

Xl ' 3 + |hi-i,j| 

(3.24) 

_ \{hi,j+i ~ 2/ijj + /ijj-i)| 
v '' 3 \hij+i\ + + \hij-i\ 

At points where h it j- 1 does not exist, 

(3.25) 

_ K^y+i ~ h*j)\ 

V '' 3 

(3.26) 

and where /ijj + 1 does not exist, 


\{hi,j-l ~ hi,j)\ 
y '' } -f \hij\ 

(3.27) 


The scaling parameter v and the dissipation constant k (the value of which is 
assigned empirically) are used to establish the level of artificial dissipation to be 
introduced. 



= k, max(i/ Si _ l j , v Xij ) 


(3.28) 
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= K max K SlJ -i . "vt, ,) (3-29) 

The corrected values of any time-stepped variable /(= h,uh,etc.) at the new 
time step are then given by the following equation : 



mod 


= f k+l 


+ 


- /«') - - f,%) l 

+K, +i Wu+1 - /«') - (/« ' 1 - /Sj + - t)l (3.30) 


where subscript mod refers to the modified value of /. The dissipation constant /c 
is used to regulate the amount of dissipation. This procedure is equivalent to adding 
second-order dissipative terms to the original governing equations. The actual nu- 
merical eddy viscosity coefficient is of the order of ku x Ax 2 / At (in ^-direction). This 
indicates that the influence of k on the results depends upon the gradients in the 
flow depth as well as grid size. The influence of k in smooth regions is minimal 
since v tends to zero in such cases. The value of k is chosen such that it is as small 
as possible but at the same time smoothens the high-frequency oscillations. The 
presence of the friction term , which is dissipative in nature, tends to inhibit the 
oscillations. Therefore, a smaller value of k is taken in those situations. 


3.6 Weighted Upwind-Downwind Method 


In the proposed method, first-order upwinding is used in the predictor part while 
first-order downwinding is used in the corrector part. Also, equal weightages are 

{ 

assigned to both predicted and corrected values while determining the values of the j 
variables at the new time level (see Eq. 3.16). Assigning of equal weightages to : 
upwinding and downwinding results in second-order accuracy and the consequent 
dispersive errors. Therefore, in the present study, it is proposed to determine the j 
values at the time level t + At as 1 



^*, + ( 2 - 0 ,)^; 

2.0 


(3.31) j 
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where w is a weighting coefficient whose value is greater than one but less two. 
u= 1 results in the second-order accurate scheme while oj = 2 results in a first-order 
upwind scheme without dispersive errors. In such a case, smoothening of the solution 
by Jameson’s method as described in the previous section may not be necessary. If 
we choose 1 < u> < 2 , we can optimise the value of u > , so that we minimize the loss 
of accuracy while at the same time removing dispersive errors. 


3.7 closure 

In this chapter, details of the Overlapping Control Volume (OCV) method for the 
solution of the shallow water flow equations are presented. The implementation 
of the boundary conditions using the reflection technique is discussed. A weightage 
upwind-downwind method is proposed to obtain less than second-order accurate so- 
lutions, but which have minimum dispersive errors. The proposed scheme is verified 
in the next chapter. 
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Chapter 4 


Verification of the Model 


The numerical technique presented in chapter 3 is verified here by applying the 
technique to obtain the flow profiles in channel transitions. The numerical results 
are compared with the corresponding results obtained by Bhallamudi and Chaudhry 
(1992) with a finite difference method. The numerical results are also compared with 
the laboratory test data available in the literature. Results for the two test cases 
are discussed in the following sections. 


4.1 Symmetrical Straight Wall Contraction 

Coles and Shintaku(Ippen et al. 1951) tested supercritical flow in the contraction 
shown in Fig.4.1. The upstream depth, ho , was 0.0305m and the upstream Froude 
number, Fo , was equal to 4.0. These test results were simulated in the present 

study using a grid Ax = 0.02415m and 12 grid divisions in the y-direction. The j 

| 

dissipation constant k =0.8 and the Courant number(C^) was equal to 0.80. It was j 
assumed that the friction and bottom slopes were equal to zero. A uniform depth of I 
0.0305m, streamwise velocity of 2.188 m/s and zero transverse velocity were specified j 

at every grid point as the initial conditions. During the computations, h, u and v j 

I 

at the upstream section were specified as 0.0305m , 2.188m/s and zero, respectively, j 

j 

to simulate the upstream boundary conditions. No condition was specified at the ; 
downstream boundary; however, values of the variables at the downstream end were j 
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extrapolated from the interior points. The flow conditions were computed upto a 
time when the flow became steady. 

Fig.4.2 and Fig.4.3 show the comparison between the results obtained using 
the OCV technique and the results obtained by employing finite difference method 
(Bhallamudi and Chaudhry, 1992). Also shown in these figures are the measured 
values obtained by Coles and Shintaku (Ippen et al. 1951). Fig. 4.2 shows the 
water surface profile along the centerline while Fig. 4.3 shows the water surface pro- 
file along the wall. It is evident from these figures that the agreement between the 
results obtained using the OCV technique and the results obtained using the finite 
difference method (Bhallamudi and Chaudhry, 1992) is very good. Therefore, it can 
be concluded that the OCV technique is as good as the finite-difference method for 
solving the shallow water flow equations. The OCV technique, moreover, does not 
require a coordinate transformation when applied to non-rectangular physical do- 
mains. This is the main advantage of the OCV technique. The agreement between 
the experimental and computed water surface profiles is good along the walls and 
at the centerline where the flows are smooth. This is not the case for the centerline 
profile in the vicinity of strong shocks. However, the computed results may be con- 
fidently used for engineering applications e.g. for selecting wall height (Bhallamudi 
and Chaudhry, 1992) . 

To investigate the effect of artificial viscosity on the computed results, the OCV 
technique is applied to the present case with k values of 0.4 and 0.8 . The computed 
water surface profiles along the centerline for these two runs are shown in Fig 4.4 
. It is clear from this figure that a two-fold increase in the artificial viscosity term 
changes the results only slightly in the vicinity of the shock, but does not affect the 
results at other sections. - 

To investigate the effect of grid size on the computed results, the OCV technique 

_ j 

is applied to the present case by taking a grid size which is two times the grid size j 
taken earlier. The dissipation constant k for this run is equal to 0.8 . The results ; 
along the centerline for this case are shown in Fig. 4.5 . It can be seen from Fig. 

4.5 that a two-fold increase in the grid size smears the shock and also reduces the j 
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maximum height slightly. As noted earlier in Chapter-3 , the amount of dissipation 
introduced by the artificial viscosity procedure increases with an increase in grid 
size, and hence the smearing of shock. 


4.2 Super critical flow in a gradual expansion 


Rouse et al.(1951) suggested the following equation for the boundaries of a gradual 
expansion (Fig 4.6 ) 

if* = (tAt) 1 ' 5 + 1 (4.1) 


A transition with the above boundaries is gradual enough so that it reduces 
the short wave effects and prevents flow seperation. The OCV technique was used 
to simulate supercritical flow in an expansion with the upstream Froude number, 
Fq= 2.0 and the upstream depth to width ratio, ho/bo= 0.25. Steady supercritical 
flow conditions at the upstream section were specified at every grid point as the 
initial condition. The values of h, u and v specified at upstream section as the 
boundary conditions were 0.0305m, 2.288m/s and zero respectively. Manning’s n 
was specified as 0.012 and the bottom slope was zero. Computations were performed 
with a Courant number equal to 0.90, dissipation constant «=0.3 , and a grid size of, 
Ax=0. 01525m. Seven grid divisions are taken in y- direction . Fig. 4.7 and Fig. 4.8 


show the comparison between the results obtained by the OCV method and those 
obtained by employing the finite difference method (Bhallamudi and Chaudhry, 
1992) . Fig 4.7 shows the water surface profile along the centerline while Fig. 
4.8 shows the results along the wall. Experimental data obtained by Rouse et al. 
(1951) are also shown in these figures. It can be .seen from these figures that the 
agreement between the results obtained using the OCV technique and the results 
obtained using the finite-difference method (Bhallamudi and Chaudhry, 1992) is 
good. The numerical results also match with the experimental data (Rouse et al. 
1951) satisfactorily. These results confirm that the OCV technique is as good as the 
finite-difference method for solving the shallow water flow equations. The numerical 
results matched satisfactorily with the experimental results because the assumption 
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of hydrostatic pressure distribution is valid in this case. It can be seen that, unlike 
in a contraction, the water surface profile in Rouse’s expansion is very smooth. 

The OCV technique with weighted upwinding-downwinding (without the use of 
artificial viscosity procedure) is also used to simulate the experiments in the Rouse’s 
expansion. The computed water surface profiles along the centerline and along the 
wall for this case are presented in Figs. 4.9 and 4.10 , respectively. The value of 
the weightage parameter, to was equal to 1.31 in these runs. It is evident from 
these figures that the weighted upwind-downwind method gives as good results as 
the finite-difference method. However, it should be mentioned here that the results 
obtained using the weighted upwind-downwind method are very sensitive to the 
value of to. Large values of to results in excessive numerical dissipation. The to 
value was chosen such that it is as close to unity as possible while at the same 
time eliminating higher order numerical oscillations. It should also be noted here 
that satisfactory results could not be obtained for the case of contraction when 
the weighted upwind-downwind method was used. Further studies are needed to 
improve the accuracy of this technique. 

The effect of k on the computed results for the expansion was investigated by 
taking two different values of k equal to 0.3 and 0.6 . The results for the centerline 
and the wall are presented in Figs. 4.11 and 4.12 , respectively. It can be ob- 
served that a two-fold increase in the artificial viscosity does not affect the results 
significantly. 
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Fig 4.2 Water surface profiles along the 
Centerline of Contraction of Fig 4.1 
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Fig 4.3 Water surface profile along the wall 
for Contraction of Fig 4.1 
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Fig 4.4 Effect of artificial viscosity 
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Fig 4.7 Water surface profile along the 
Centerline for expansion of fig 4.6 
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Fig 4.8 Water surface profile along the 
Wall for Expansion of Fig 4.6 
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Fig 4.9 Water surface profile along the centerline 
(Weighted Upwind-Downwind Method) 
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Fig 4.10 Water surface profile along the Wall 
(Weighted Upwind-Downwind Method) 
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Fig 4.1 1 Effect of artificial viscosity for expansion 
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problem (Centerline profile ) 
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Fig 4.12 Effect of artificial viscosity for 


expansion problem (Wall profile) 




Chapter 5 


Summary and Conclusions 


In this study, a control volume technique has been developed for solving the shallow 
water flow equations which describe the unsteady open-channel flows. The proposed 
control volume technique is based on the principle of overlapping control volumes 
suggested by Verma and Eswaran (1996) . The time discretization is done us- 
ing a MacCormack type Predictor-Corrector approach. First-order upwinding and 
first-order downwinding are used for computing fluxes at the control surface during 
predictor and corrector calculations, respectively. The higher-order dispersive er- 
rors which result because of the second— order accuracy of the scheme are eliminated 
using the artificial viscosity approach. Alternatively, a weighted upwind-downwind 
method is proposed for this purpose. 

The OCV technique has been applied for simulating the steady two-dimensional 
supercritical flows in channel transitions. The unsteady computations were used 
in a false-transient approach to obtain the steady state solutions. Both channel 
expansions and contractions were considered as test cases. The numerical results 
obtained using the OCV technique matched very well with the earlier numerical 
solutions obtained using the finite-difference approach. The numerical results also 
matched satisfactorily with the experimental data for those cases where the shallow 

water flow assumption is valid. 

The main advantage of the proposed OCV technique is that it can be applied to 
general non-orthogonal domains without having to do the cumbersome coordinate 
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transformation. Unlike the finite-element techniques, it is easy to implement and 
the computational requirements are are not high. As compared to the usual finite- 
volume techniques, the proposed OCV technique uses the grid point coordinates, 
directly, to form the control volumes. It does not require the determination of the 
intermediate points for forming the control volumes. 

The proposed weighted upwind-downwind method worked well for the case of 
channel expansions in which the water surface profile was smooth. However, it did 
not give correct results for the case of channel contractions in which shocks are 
present . OCV technique with artificial viscosity approach worked well for both the 
cases. However, artificial viscosity approach has been controversial in recent years 
since it uses an adjustable parameter. Further studies- are required to incorporate 
the Total Variation Diminishing (TVD) principle into the OCV technique. 



REFERENCES 


Akanbi, A. A. and Katapodes, N.D. (1988) “Model for flood propagation on 
initially dry land” , Jour, of Amer. Soc. Civ. Engrs. , Vol.114, (7), pp689-706. 

Alcrudo, F. Garcia-Novvarro, P. and Saviron, J.M. (1992) “Flux-difference 
splitting for 1-Dimensional open channel flow equations” . Inti. Jour, for numerical 
methods in fluids, Vol 14, ppl009-1018. 

Alcrudo, F. and Garcia, N.P. (1993) A TVD scheme in Finite— volumes for sim- 
ulation of two-dimensional discontinuous flows. Inti, Jour, for Numerical Methods in 
Fluids. Vol. 16, pp489-505. 

Anderson, D.A., Tannehill, J.D. and Pletcher, R.H., (1984), Computational 
Fluid Mechanics and Heat transfer , McGraw-Hill, New York. 

Bellos, C.V. , Soulis, J.V. and Sakkas, J.G. (1991) Computation of two- 
dimensional clambreak induced flows. Computational Mechanics Publications. Advanced 
Water Res. Vol. 14, No. 1 . 

Bhallamudi, S.M. and Chaudhary, M.H., (1992), ’’Computation of flows in. 
Open-Channel Transitions”, Journal of Hydraulic Research, IAHR, Vol. 30, No.l, pp. 
77-93. 

Casulli, V. (1990), “Numerical Simulation ofehallow water flow ”, in Computational 
methods in Surface Hydrology, (Ed: Gambolti, G. et al.), Springer^-Verlag. 

Chaudhary, M.H. (1993), Open-Channel Flow, Prentice hall, New Jersey. 

Chaudhary, M.H. , Bhallamudi, S.M. and Gharangik, A.M. (1985) em Math- 
ematical modelling of flows in Baker Bay - Technical report prepared for U.S.Army 
Corpse of Engineers, Portland, Oregon. 

Dammular, D. , Bhallamudi, S.M. and Chaudhry, M.H. (1989) Modelling of 
Unsteady flow in Curved channels, Jour, of Hydr. Engg., ASCE, 115 (11) ; 1479-95. 

Fennema, R.J. and Chaudhry, M.H. (1989) Implicit methods for two-dimensional 
unsteady free surface flows, Jour, of Hydr. Res. , IAHR, 27(3) ; pp 321-32. 

Fennema, R.J. and Chaudhry, M.H. (1990) Numerical Solution of Two-Dimensional 
Transient Free Surface Flows, Jour, of Hydr. Eng., Amer. Soc. Civ. Engrs., Vol 116, 
Aug., ppl013-1034. 

Galland, J.C. et al. (1991) TELEMAO, A New Numerical Model for Solving 
Shallow Water Flow Equations. Adv. Water Res. Vol. 14, No.3 . 



MacCor mack, R.W., (1969). The Effect of viscosity in hypervelocity impact era 1 
tering, Amer. Inst. Aeronautics and Astronautics, Paper 69-354, Cincinnati, Ohio. 

Nujic, M. (1995) Efficient implementation of non-oscillatory schemes for the com- 
putation of free-surface flows, IAHR, Vol.33, No.l, ppl01-110. 

Peyret, R. and Taylor, T.D. (1983) Computational Methods for Fluid Flow, 
Springer-Verlag, New York, 1983. 

Puri, A.N. and Kao, C.Y., (1985), “Numerical Modelling of Subcritical Open 
Channel Flow Using The K-e Turbulence Model and The Penalty Function Finite Element 
Technique,” Appl. Math. Modelling, vol. 9, No. 2, Apr., pp. 82-88. 

Putti, M , Yell, W.W.G. and Mulder, W.A. (1990) ‘A Tringular finite volume 
approach with high resolution upwind terms for the solution of ground water transport 
equation’ , Water Res. Research, 26, pp2865-2880. 

Roache, P.J., (1972), Computational Fluid Dynamics, Hermosa Publishers. 

Rouse, H., Bhoota, B.V., and Hsu, E.Y. (1951). Design of channel espansions, 
Trans., Amer. Soc. Civil Engrs., 116:347 . 

Savic, L.J. and Holly Jr., F.M. (1993) Dambreak flood waves computed by mod- 
ified Godunov method. , IAHR, Vol. 31 , No.2 , ppl87-204. 

U.S. Army Corpse of Engrs. (1985) “TABS-2 : User’s Manual for the Gener- 
alized Computer program system, open channel flow and sedimentation” , Waterways 
Experiment Station, Vicksburg. 

Verma, A.K , Eswaran, V. (1996) ‘Overlapping Control Volume approach for 
convection-diffusion equation’, Int. J. Numer. Methods in Fluids [in press]. 

Verma, A.K. , Bhallamudi, S.M. and Eswaran, V. (1996) ‘Overlapping Control 

Volume Method for Solute Transport’ Paper Communicated to Inti. Journal for Numer. 

Methods in Engg. . 

Yang, J.Y. , Hsu, C.A. and Chang, S.H. (1993) Computation of free surface 
flows, Part 1 : One-dimensional dambreak flow, IAHR, Vol.31, ppl9-34. 



